library(readr)
library(dplyr)
library(ggplot2)
library(lubridate)
library(tidyr)
library(DescTools)
library(gridExtra)

1 Analisis preliminar de intensidad

Cargamos una base generada por la unión de todas las bases mensuales únicamente quedandonos con las variables de interes:

datos <- read_csv("/Users/pablogandia/Desktop/trafico_meses_juntos.csv")
head(datos)

Vemos que esta base de datos contiene todas las zonas con todos los meses, con granularidad de 5 minutos. Ha tardado aproximadamente 20 minutos en cargar y ocupa 9GB de RAM:

unique(datos$Mes)

##  [1]  1  2  3  4  5  6  7  8  9 10 11 12

1.1 Ejemplo de intensidad en la carretera del Politécnico (Av. Naranjos)

datos <- datos %>%
  mutate(
    Mes = month(FechaHora),
    Dia = day(FechaHora)
  )

arc <- datos %>%
  filter(IdPM == 2046, Mes == 1, Dia > 1, Dia < 8)

ggplot(arc, aes(x = FechaHora, y = Intensidad)) +
  geom_line(linewidth = 1.2, color = "steelblue") +
  geom_point(size = 1.5, color = "black") +
  labs(
    title = "Intensidad durante la primera semana en zona 2046",
    x = "Fecha y hora",
    y = "Intensidad"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 14)
  )

Vemos que no tiene sentido que cada 5 minutos tome valores tan extremandamente altos, consultamos con la empresa proveedora y nos confirma que la intensidad son aproximiaciones de la media diaria, es decir, para conocer exactamente la cantidad de coches que pasan en esos 5 minutos deberíamos dividir intensidad entre 12 (60/5).

datos <- datos %>%
  mutate(IntensidadReal5mins = (Intensidad / 60) * 5)

ordenado <- datos %>%
  arrange(desc(IntensidadReal5mins))

head(ordenado,5)

Observamos zonas con Intensidades reales cada 5 minutos completamente desorbitadas, por lo que planteamos que pueden darse fallos en las espiras. Observemos estos datos más de cerca:

arc <- datos %>%
  filter(IdPM == 1418, Mes == 8, Dia == 2)

ggplot(arc, aes(x = FechaHora, y = IntensidadReal5mins)) +
  geom_line(linewidth = 1.2, color = "steelblue") +
  geom_point(size = 1.5, color = "black") +
  labs(
    title = "Intensidad rara zona 1418",
    x = "Fecha y hora",
    y = "Intensidad (5 min)"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 14)
  )

Efectivamente, de repente se dan picos muy diferenciados de sus datos vecinos, y eso que son intervalos de 5 minutos donde es imposible que pasen tantos coches.

datos <- datos %>%
  mutate(
    Mes = month(FechaHora),
    Dia = day(FechaHora)
  )

arc <- datos %>%
  filter(IdPM == 1612, Mes == 3, Dia == 6)

ggplot(arc, aes(x = FechaHora, y = IntensidadReal5mins)) +
  geom_line(linewidth = 1.2, color = "steelblue") +
  geom_point(size = 1.5, color = "black") +
  labs(
    title = "Intensidad rara zona 1612",
    x = "Fecha y hora",
    y = "Intensidad (5 min)"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 14)
  )

1.2 Decisión de función agregada

Como nuestra intención es pasar la granularidad a datos horarios, hay que encontrar la medida agregada que mas respete los datos, pero que a la vez nos proteja frente a estos grandes picos, que claramente son errores ya que son físicamente imposibles de darse.

Paro esto comparamos medidas agregadas, como la media, mediana y media winsorizada al 5 y al 10:

wmean_func <- function(x, p = 0.05) {
  x <- x[!is.na(x)]
  q_low <- quantile(x, p)
  q_high <- quantile(x, 1 - p)
  x[x < q_low] <- q_low
  x[x > q_high] <- q_high
  mean(x)
}

hourly <- datos %>%
  group_by(IdPM, Hora, Mes, Dia) %>%
  summarise(
    median = median(Intensidad, na.rm = TRUE),
    mean   = mean(Intensidad, na.rm = TRUE),
    wmean_05pct = wmean_func(Intensidad, p = 0.05),
    wmean_10pct = wmean_func(Intensidad, p = 0.10),
  )

datos <- datos %>%
  left_join(hourly, by = c("IdPM", "Hora", "Mes", "Dia"))

head(datos)

Volvamos a ver el comportamiento con los datos que antes habian sido extraños

arc <- datos %>%
  filter(IdPM == 1612, Mes == 3, Dia == 6)

arc_largo <- arc %>%
  select(FechaHora, Intensidad, median, mean, wmean_05pct, wmean_10pct) %>%
  pivot_longer(
    cols = -FechaHora,
    names_to = "Tipo",
    values_to = "Valor"
  )

ggplot(arc_largo, aes(x = FechaHora, y = Valor, color = Tipo)) +
  geom_line(linewidth = 1.2) +
  geom_point(data = arc_largo %>% filter(Tipo == "Intensidad"),
             aes(x = FechaHora, y = Valor), color = "darkblue", size = 1.5) +
  labs(
    title = "Intensidad 5 minutos a horaria",
    x = "Fecha y hora",
    y = "Intensidad",
    color = "Tipo"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 14)
  )

Vemos que todas las otras medidas se dispersan mucho, pero la mediana en cambio sigue la distribución acorde a sus datos anteriores, por lo que parece una función agregada perfecta para cambiar la granularidad.

Veamos como funciona en una distribución sin valores erroneo, como puede ser de nuevo la avenida de enfrente del politécnico:

arc <- datos %>%
  filter(IdPM == 2046, Mes == 1, Dia == 1)


arc_largo <- arc %>%
  select(FechaHora, Intensidad, median, mean, wmean_05pct, wmean_10pct) %>%
  pivot_longer(cols = -FechaHora, names_to = "Tipo", values_to = "Valor") %>%
  mutate(Tipo = recode(Tipo,
                       Intensidad     = "Original",
                       median         = "Mediana",
                       mean           = "Media",
                       wmean_05pct    = "Winsor 5%",
                       wmean_10pct    = "Winsor 10%"))

ggplot(arc_largo, aes(x = FechaHora, y = Valor, color = Tipo)) +
  geom_line(linewidth = 1.2) +
  geom_point(data = arc_largo %>% filter(Tipo == "Original"),
             aes(x = FechaHora, y = Valor), color = "darkblue", size = 1.5) +
  labs(
    title = "Intensidad 5 minutos a horaria",
    x = "Fecha y hora",
    y = "Intensidad",
    color = "Serie"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 14)
  )

Vemos que restan un poco de variabilidad, pero en general funcionan muy similares, la media y media winsorizada al 10 si que despuntan ligeramente en los picos de intensidad, cuyo dato puede ser correcto, pero la debemos sacrificar cierta variabilidad para poder cambiar de granularidad, a la vez que evitamos los errores momentaneos.

Veamos como tratan los atípcios en general, para esto seleccionamos 10 zonas aleatorias, ya que si no sería imposible interpretar los atípicos, ya que como hemos visto en otros análisis, al ser las zonas tan diferentes, y variar tanto el tráfico dependiendo de la noche y el dia, hay una cantidad enorme de valores atípicos:

idpm_seleccionados <- sample(unique(datos$IdPM), 10)
data2 <- datos %>%
  filter(IdPM %in% idpm_seleccionados)

data2

variables <- c("median", "mean", "wmean_10pct", "Intensidad")

for (var in variables) {
  
  df_plot <- data2 %>%
    select(Hora, !!sym(var)) %>%
    filter(!is.na(Hora), !is.na(!!sym(var))) %>%
    rename(valor = !!sym(var))
  
  p <- ggplot(df_plot, aes(x = factor(Hora), y = valor)) +
    geom_boxplot(outlier.shape = 1, fill = "skyblue", color = "black") +
    labs(
      title = paste("Boxplot de", var, "por hora"),
      x = "Hora del día",
      y = var
    ) +
    theme_minimal()
  
  print(p)
}

Vemos que aunque si bien las otras medidas de agrupación no contienen valores tan imposibles como los que se dan en la intensidad cada 5 minutos, la mediana es la medida más robusta de agrupación para la intensidd,y por tanto la que usaremos para pasar a datos horarios.

Veamos si varia mucho el comportamiento general generando la simulacion de lo que sería la agrupación horaria con la mediana:

cols_to_drop <- c(
  "Ocupacion", "Velocidad",
  "FiabilidadIntensidad", "FiabilidadOcupacion", "FiabilidadVelocidad",
  "Minuto",
  "DiaSemana"
)

df_reducido <- datos %>%
  select(-any_of(cols_to_drop))

df_horario <- df_reducido %>%
  group_by(IdPM, Año, Mes, Dia, Hora, FechaHora) %>%
  summarise(Intensidad_mediana = median(Intensidad, na.rm = TRUE), .groups = "drop") %>%
  arrange(FechaHora) %>%
  select(FechaHora, IdPM, Intensidad_mediana)

head(df_horario)

p1 <- ggplot(data2, aes(x = Intensidad)) +
  geom_histogram(bins = 50, fill = "steelblue", color = "white") +
  labs(
    title = "Distribución de Intensidad (5 mins)",
    x = "Intensidad (vehículos/hora)",
    y = NULL
  ) +
  theme_minimal() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

p2 <- ggplot(df_horario, aes(x = Intensidad_mediana)) +
  geom_histogram(bins = 50, fill = "steelblue", color = "white") +
  labs(
    title = "Distribución de Intensidad (mediana horaria)",
    x = "Intensidad (vehículos/hora)",
    y = NULL
  ) +
  theme_minimal() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

grid.arrange(p1, p2, ncol = 2)

media1 <- datos %>%
  group_by(Hora) %>%
  summarise(Intensidad = mean(Intensidad, na.rm = TRUE))

df_horario$Hora <- hour(df_horario$FechaHora)

media2 <- df_horario %>%
  group_by(Hora) %>%
  summarise(Intensidad_mediana = mean(Intensidad_mediana, na.rm = TRUE))

p1 <- ggplot(media1, aes(x = Hora, y = Intensidad)) +
  geom_line(color = "steelblue", linewidth = 1) +
  geom_point(color = "steelblue", size = 2) +
  labs(title = "Intensidad por hora") +
  theme_minimal()

p2 <- ggplot(media2, aes(x = Hora, y = Intensidad_mediana)) +
  geom_line(color = "steelblue", linewidth = 1) +
  geom_point(color = "steelblue", size = 2) +
  labs(title = "Intensidad mediana por hora") +
  theme_minimal()

grid.arrange(p1, p2, ncol = 2)

Concluimos con que la mediana si es una buena medida robusta y representativa para hacer la agrupación horaria, al menos en la muestra seleccionada.

2 Ocupación preliminar

Ahora buscaremos una función óptima para la variable ocupación

arc <- datos %>%
  filter(
    IdPM == 2046,
    Mes == 1,
    Dia == 1,
  )

ggplot(arc, aes(x = FechaHora, y = Ocupacion)) +
  geom_line(linewidth = 1.2, color = "steelblue") +
  geom_point(size = 1.5, color = "black") +
  labs(
    title = "Ocupación durante un dia en zona 2046",
    x = "Fecha y hora",
    y = "Ocupación"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    plot.title = element_text(size = 14)
  )

Vemos que aunque presenta algunos picos, la gran mayoria de las veces desciende y fluctua mucho, por eso planteamos que quizas una mediana no fuera tan buena, porque estaríamos sacrificando la variabilidad de esta variable tan cambiante que funciona por picos.

wmean_func <- function(x, p = 0.05) {
  x <- x[!is.na(x)]
  q_low <- quantile(x, p)
  q_high <- quantile(x, 1 - p)
  x[x < q_low] <- q_low
  x[x > q_high] <- q_high
  mean(x)
}

hourly <- datos %>%
  group_by(IdPM, Hora, Mes, Dia) %>%
  summarise(
    median_OC = median(Ocupacion, na.rm = TRUE),
    mean_OC   = mean(Ocupacion, na.rm = TRUE),
    wmean_05pct_OC = wmean_func(Ocupacion, p = 0.05),
    wmean_10pct_OC = wmean_func(Ocupacion, p = 0.10),
  )

datos2 <- datos %>%
  left_join(hourly, by = c("IdPM", "Hora", "Mes", "Dia"))

head(datos2)

variables <- c("median_OC", "mean_OC", "wmean_10pct_OC", "Ocupacion")

for (var in variables) {
  
  df_plot <- datos2 %>%
    select(Hora, !!sym(var)) %>%
    filter(!is.na(Hora), !is.na(!!sym(var))) %>%
    rename(valor = !!sym(var))
  
  p <- ggplot(df_plot, aes(x = factor(Hora), y = valor)) +
    geom_boxplot(outlier.shape = 1, fill = "skyblue", color = "black") +
    labs(
      title = paste("Boxplot de", var, "por hora"),
      x = "Hora del día",
      y = var
    ) +
    theme_minimal()
  
  print(p)
}

3 Velocidad

Puesto que como nos indicaron, se utiliza la media, tambien la utilizaremos para encontrar el computo horario.

Y esque los valores de velocidad, sin contar los errores, están limitados a la características de la vía: Por mucho que se busque, las carreteras están limitadas a un tope de velocidad.

4 Generar las variables horarias:

Para respetar la mayor parte de la variabilidad decidimos seleccionar las siguientes variables:

Intensidad: Mediana, Máximo y Desviación Típica Ocupación: Media y Máximo Velocidad: Media Fiabilidades: Si durante la hora ha habido un momento de fiabilidad menor de 0 se establece FalloVariable como True